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Abstract 

We study the spectral gap of the Wilson-Dirac operator in two-flavour lattice QCD as a 
function of the lattice spacing a , the space-time volume V and the current-quark mass to. 
It turns out that the median of the probability distribution of the gap scales proportionally 
to to and that its width is practically equal to a/\/V . In particular, numerical simulations 
are safe from accidental zero modes in the large-volume regime of QCD. 


1. Introduction 

The formulation of lattice QCD proposed by Wilson long ago preserves many sym¬ 
metries of the continuum theory exactly [1]. An infamous exception are the chiral 
symmetries, and although the symmetry-violating terms vanish proportionally to 
the lattice spacing (or its square if the theory is Symanzik-improved [2,3]), the pres¬ 
ence of these terms complicates the lattice theory considerably. In particular, the 
fact that the massive Wilson-Dirac operator is not protected from arbitrarily small 
eigenvalues may lead to instabilities in numerical simulations. 

Such instabilities were not observed, however, in recent simulations of the two- 
flavour Wilson theory on large lattices, not even at the smallest quark masses consid- 
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ered [4-6]. It is tempting to attribute the absence of instabilities in these simulations 
to the use of a new simulation algorithm, but this explanation cannot be right, be¬ 
cause the distribution of the spectral gap of the lattice Dirac operator (and thus 
the probability to find exceptionally small eigenvalues) is determined by the lattice 
action and the functional integral, and not by the method used to evaluate the latter. 

Our aim in this paper is to clarify the situation by calculating the distribution of 
the gap on a set of lattices, using numerical simulations. In particular, we wish to 
determine, as explicitly as possible, the range in parameter space, where the Wilson 
theory can be simulated without running into instabilities. 


2. Algorithmic stability and the spectral gap 

In this section we introduce our notation and discuss the relevance of the distribution 
of the spectral gap for the stability of lattice QCD simulations. For any unexplained 
notation see ref. [3]. 

2.1 Lattice Dirac operator 

The lattice theory is set up as usual on hypercubic lattices with spacing a. Periodic 
boundary conditions are imposed on all fields and in all directions, except for the 
quark fields, which are taken to be antiperiodic in time. Throughout this paper, we 
assume that there is a doublet of sea quarks with equal mass, although many results 
are likely to remain valid if the number of quark flavours is larger than two. 

While the action of the gauge field will always be the Wilson plaquette action, we 
shall consider various lattice Dirac operators D m (the subscript indicates that D m 
includes the quark mass term). The Wilson-Dirac operator 


D m — H - m 0 , 

(2.1) 

D w = \ { 7/1 (V* + V„) - aV*V M } , 

(2.2) 


is the one we are primarily interested in, as well as its 0(a)-improved version [2,3]. 
In these equations, and denote the gauge-covariant forward and backward 
difference operators, and mo the bare quark mass. 

Rather than the Dirac operator itself, we prefer to consider the hermitian operator 

Qm = 75-tkn (2-3) 
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in the following. The determinants of these two operators are the same, but the fact 
that the spectrum of Q m is real simplifies the discussion considerably. On a finite 
lattice, and for any specified gauge field, we then define the spectral gap 

p = min {|A| | A is an eigenvalue of Q m } • (2-4) 

Evidently, the gap is a well-defined function of the gauge field and so is the spectral 
asymmetry 

r 1 = \{N + -N_}, (2.5) 


where N± are the numbers of positive and negative eigenvalues of Q m . 

In formulations of lattice QCD that preserve chiral symmetry via the Ginsparg- 
Wilson relation [7-11], the gap is bounded from below by the bare current-quark 
mass m and the asymmetry vanishes if m > 0 (see ref. [12], for example). However, 
these properties are not guaranteed in the Wilson theory, and it is possible, in the 
presence of some gauge-field configurations, that the gap is much smaller than m 
and that the asymmetry assumes a non-zero value. 

2.2 Sources of instability 

In the large-volume regime of QCD, and at large quark masses, the probability dis¬ 
tribution p{n) of the gap typically looks like the one shown in the upper plot in fig. 1. 
The subsets of gauge fields, where p is far below the central value of the distribu¬ 
tion, occur with such a small probability in this case that their contributions to the 
common physical observables can be safely neglected. Numerical simulations, using 
a preconditioned Hybrid Monte Carlo (HMC) algorithm [13], for example, will then 
normally run smoothly and produce a representative sample of field configurations 
as expected. 

When the quark mass decreases, the centre of the gap distribution moves towards 
smaller values and eventually is no more than one or two standard deviations away 
from the origin (lower plot in fig. 1). At this point numerical simulations may run 
into instabilities for the following reasons: 

(1) Integration instabilities. The HMC algorithm obtains the next field configu¬ 
ration by integrating the appropriate molecular-dynamics equations in field space, 
followed by an accept-reject step at the end of the integration. It is possible that 
the molecular-dynamics trajectories pass through field configurations where the gap 
p is exceptionally small. The numerical integration becomes unstable in this case 
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Fig. 1. Qualitative form of the normalized probability distribution p(p) of the gap 
(full lines), at two values of the quark mass (upper and lower plot). The normalized 
weighted distributions proportional to p(p)/p 2 are also shown (dotted lines). 

and liable to rounding errors. In particular, the reversibility of the integration is 
then no longer guaranteed, which invalidates the algorithm. 

(2) Ergodicity. Along the molecular-dynamics trajectories, the spectral asymmetry 
rj tends to be preserved. Changes in the asymmetry are in fact excluded in the limit 
where the integration step-size goes to zero. However, the sectors in field space, 
in which the asymmetry does not vanish, may become statistically relevant in the 
situation considered here. The HMC algorithm is then likely to give wrong results, 
because it may get stuck, for a very long time, in any one of these sectors. 

(3) Sampling inefficiencies. Observables that are sensitive to the small eigenvalues 
of the Dirac operator, such as the pion propagator, may be poorly sampled by the 
simulation. The effect is illustrated by the dotted curves in fig. 1, which represent the 
distributions of the gap reweighted by the “observable” 1 //x 2 . Expectation values of 
such quantities are difficult to compute reliably, since the subspace of fields where 
p is very small is only rarely visited in the course of the simulation. 

Instabilities of this kind can lead to underestimated statistical errors and incorrect 
results, and they may even suggest the presence of a phase transition when there is 
none. Since the QCD functional integral remains well defined, also in this difficult 
regime, improved simulation techniques may conceivably be developed, which do not 
suffer from any instabilities. Whether this is worth the effort is not obvious, however, 
because the theory may be strongly affected by lattice artefacts (resulting from a 
competition of mass and discretization effects) in the critical range of parameters. 
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3. Numerical studies 


In the presence of any given gauge field, the spectral gap of the Dirac operator can be 
computed numerically, using a suitable iterative method (appendix A). Evidently, 
the histograms of the values calculated in the course of a numerical simulation 
approximate the gap distribution up to statistical errors. In this section, we report 
the results of such numerical studies and show that the data are well described by 
a few simple empirical laws. 

3.1 Simulation parameters 

The simulations listed in table 1 are part of an ongoing study of two-flavour QCD 
in the chiral regime [6]. Technically the project is based on the use of the Schwarz- 
preconditioned HMC simulation algorithm introduced in ref. [4], which allows the 
theory to be simulated in a range of lattices and quark masses that was, in practice, 
inaccessible so far. 

Except for the last run in table 1, the unimproved Wilson theory was simulated, 
at inverse bare coupling /3 and hopping parameter k = (8 + 2amo)“ 1 . In run D\ the 
coefficient c sw of the Sheikholeslami-Wohlert improvement term [2,3] was set to the 
value determined by the ALPHA collaboration [14], A number IV c fg of statistically 
decorrelated gauge-field configurations was generated in each case and later used for 
the calculation of the gap distributions. 


Table 1. Simulation runs included in this study 


Run 

Lattice 

P 

£-SW 

K 

Acfg 

am n 

Ai 

32 x 24 3 

5.6 

0 

0.15750 

64 

0.2744(21) 

A-2 




0.15800 

109 

0.1969(16) 

A3 




0.15825 

100 

0.1554(31) 

a 4 




0.15835 

100 

0.1204(44) 


64 x 32 s 

5.8 

0 

0.15410 

100 

0.1965(8) 

b 2 




0.15440 

101 

0.1481(11) 


64 x 24 3 

5.6 

0 

0.15800 

116 

0.1986(10) 

D\ 

48 x 24 3 

5.3 

1.90952 

0.13550 

104 

0.3265(11) 
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Fig. 2. Normalized histograms of the spectral gap g, as obtained in the runs A 1 —A 4 
and Ci. The bin size is 1.5 MeV and the dotted vertical lines indicate the position of 
the median p of the distributions. 

Details of the simulations, the computation of the quark masses and other physical 
quantities will be given in a forthcoming publication [6]. Based on calculations of the 
Sommer reference scale ro [15], the lattice spacing in the Wilson theory is estimated 
to be about 0.08 and 0.06 frn at /3 = 5.6 and 5.8, respectively, while it is roughly 
0.09 fm in the improved theory at (3 = 5.3. In physical units the spatial sizes of the 
lattices in table 1 are thus close to 2 fm in all cases. 

As usual the conversion from lattice to physical units is ambiguous, but in this 
paper the latter serve for the purpose of illustration only. The pion masses quoted 
in the last column of table 1, for example, cover a range from 676 to about 294 MeV 
in the case of the runs A 1 -A 4 . 

3.2 Mass-dependence of the gap distribution 

The calculated distributions of the gap in the Wilson theory at (3 = 5.6 are plotted 
in fig. 2. A characteristic feature of these distributions is that their shape does not 
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Fig. 3. Median jl of the gap distribution on the 32 x 24 3 lattice (runs A 1 -A 4 ) as 
a function of the bare current-quark mass to. The straight line from the point at the 
largest mass to the origin is drawn to guide the eye. Statistical errors are negligibly 
small on the scale of this plot. 

show a strong dependence on the quark mass. Basically they are shifted to smaller 
values when the quark mass decreases. All these distributions are clearly separated 
from the origin, although the smallest mass may be quite close to where the unstable 
regime begins. 

The distributions are roughly symmetric about their median (dotted lines in fig. 2), 
which is practically also the point where they assume their maximal value. As shown 
in fig. 3, the median scales approximately linearly with the current-quark mass. A 
small but significant effect is seen at the lighter quark masses, where the median is 
pushed to slightly larger values with respect to the straight scaling curve. 

On the big lattices at fl = 5.8, the situation is essentially the same, although here 
only two simulations have been completed so far (see fig. 4). In particular, the quark 
masses in these two runs are in a range where a nearly perfect scaling of the median 
as a function of the quark mass is again observed. 

3.3 Statistical fluctuations of the gap 

An important qualitative result of our empirical studies is that the fluctuations of 
the gap become smaller when the four-dimensional volume V of the lattice increases. 
The simulation C\. for example, was performed at exactly the same coupling and 
quark mass as the simulation A 2 , but on a larger lattice. As can be seen from fig. 2, 
the width of the gap distribution obtained in run C\ is visibly reduced with respect 
to the one obtained in run A- 2 - 
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Fig. 4. Gap distributions and mass dependence of the median /i on the 64 x 32 3 
lattices (runs B i and B^). The bin size is the same as in fig. 2. 

A simple heuristic argumentation that may explain this effect goes as follows. Let 
U be a given lattice gauge field, /i the lowest eigenvalue of \Q m \ and ip the associated 
eigenvector normalized to unity. Consider a small random fluctuation U + 5U of the 
gauge field. To first order, the change 5/i in the gap is then given by 

5ji = a 4 y^^(x) f [5Q m ip\ (x). (3.1) 

X 

At each point x , there are eight contributions to this sum, which are proportional 
to the random variation 5U of the gauge field on the links attached there. Now if ip 
extends over the whole lattice, i.e. if ip it is not a localized mode, the size of these 
terms is proportional to a 3 /V. On average this implies (6/j) = 0 and ((<5/u) 2 ) oc a 2 /V. 
While this kind of reasoning is extremely superficial, the width a of the numerically 








Fig. 5. Width a of the gap distributions, given in units of a/s/V, as obtained in 
runs A 1 —D 1 . Note that the combination a\/V/a is dimensionless and can be computed 
directly from the lattice data, without intermediate conversion to physical units. The 
statistical errors were determined using the bootstrap method. 


determined gap distributions appears to scale in the suggested way (see fig. 5) f. An 
interesting special case is that of the distributions obtained in the runs A 2 and C \, 
where, as already mentioned, all parameters were the same except for the lattice 
volumes. The scaled widths a\fV / a , however, agree very well with each other. The 
approximate scaling law 


a ~ 


y/V 


(3.2) 


is actually consistent with all the data shown in fig. 5. 

Somewhat surprisingly, the spectral gap thus seems to behave thermodynamically, 
like a free-energy density for example, with probability distribution 


V 2 

p(aO «exp<j - —(/r- (/z» 


(3.3) 


The fact that a is proportional to the lattice spacing suggests, on the other hand, that 
the observed widths result from short-distance fluctuations of the gauge field. These 
fluctuations scale more slowly with the volume than those expected in lattice theories 

j" The standard deviation of /r is subject to potentially large statistical uncertainties, because the 
tails of the gap distributions are poorly sampled. We therefore define the width of the distributions 
through a = ^(v — u), where [u, v\ is the smallest range in fj ;, which contains more than 68.3% of 
the data. 
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that preserve chiral symmetry, where the eigenvalue distributions are universally 
computable, using chiral perturbation theory or random matrix theory. 


4. Spectral gap in infinite volume 


Analytical calculations of the gap distribution would evidently be very welcome at 
this point, but since chiral symmetry is violated in the Wilson theory, it is unclear 
how such calculations would proceed. The spectral density of the Dirac operator in 
infinite volume is somewhat more accessible and provides a useful reference for the 
situation on the lattices that can be simulated. 

4.1 Spectral density 

Let 01 , 02 ) • • • be the eigenvalues of Q 2 n , ordered in ascending order and counting 
multiplicities. In the following we will be interested in the spectral density 

P(«) = J im 77 ^( a “ “*)) C 4 - 1 ) 

V —MX) V * -' 

k> 1 

in infinite volume. It is, incidentally, possible to prove rigorously that the thermo¬ 
dynamic limit (4.1) exists, using a general argument based on a decomposition of 
the lattice into large blocks [16]. 

In the continuum theory, or if chiral symmetry is preserved by the lattice theory, 
the spectral density vanishes at o < m 2 . Moreover, for slightly larger values of a, 
the Banks-Casher relation [17] implies the asymptotic form 

P ( a ) = 2 —/ 2 + Q(!)) (4-2) 

a>m 2 7T V O; — 77?. 

where £ denotes the quark condensate. The spectral density thus has a well-defined 
threshold at a = m 2 in these cases (see fig. 6). 

In the following a working hypothesis is that there is a similar threshold a > 0 in 
the Wilson theory or its 0(a)-improved version (whichever is considered). We shall 
not need to know the shape of the spectral density in the vicinity of the threshold, 
but it must be guaranteed that p(a ) vanishes if a < a and that the density is non¬ 
zero at or immediately above a. It is important to understand that a is the point 
where the dense spectrum begins, while it is perfectly possible that the Wilson- Dirac 
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Fig. 6. Number of eigenvalues per bin and fm 4 at the low end of the spectral density 
in the continuum theory [cf. eq. (4.2)]. In this plot, m = 15 MeV and E = (250 MeV) 3 
was assumed. 

operator in finite volume has eigenvalues significantly smaller than \Jlk. Any part 
of the spectrum with less than 0(F) modes per energy bin is in fact suppressed by 
the factor 1/V in the limit (4.1). 

4-2 Adding valence quarks 

To be able to relate the spectral density to correlation functions of local fields, we 
will need to consider partially quenched QCD, where 2 N valence quarks are added 
to the theory [18,19]. The fermion action then becomes 

{ 2N+2 N \ 

E r {x)D m ^r{x ) + E \Dm4>k(x)\ 2 > , (4.3) 

r=l k= 1 ) 

where ip r ,ip r are the quark fields (2 sea quarks plus 2 N valence quarks) and 4 >k the 
pseudo-fermion fields that are required to cancel the valence-quark determinants. 
This defines a well-behaved euclidean field theory, whose renormalization can be 
expected to follow the usual rules. In particular, the obvious U(21V + 2) X U(1V) 
flavour symmetry and some less obvious graded symmetries restrict the possible 
counterterms to the naively expected ones. 

We now also introduce the bare scalar and pseudo-scalar densities 

S rs = VvV’s, Prs = i’r'y^s- (4.4) 

These fields renormalize in the standard manner, i.e. for unequal flavours r, s , the 
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renormalized operators are Z^S rs and ZpP rs , where the renormalization constants 
Zs and Zp are flavour-independent and can be taken to be the same as the ones in 
the theory without valence quarks. 

4-3 Resolvent & moments of the spectral density 
The resolvent 

R(z)= [°°da 2 f a) (4.5) 

Ja a z (z-a) 

is an analytic function of z with a cut along the real axis from z = a to some value 
proportional to 1/a 2 . For later convenience, a factor or 2 is included in the integral 
to improve its convergence at large eigenvalues. Evidently, the spectral density is 
proportional to the discontinuity across the cut, and so can be uniquely recovered 
from the resolvent if the latter is known [20,21]. 

For \z\ < a, the resolvent may be expanded in a convergent power series, 

°° /»00 / \ 

*) = £>***, M k = - dagg 

fc =0 Ja 

with coefficients M^ that can be written as 

Mk = a 4 ” -4 <Pi 2 (si)i , 23(*2)...i , ni(0)), n = 2k + 6. (4.7) 

Xl. •jXn — 1 

This formula assumes that there are at least n quarks, but since each moment M& 
may be considered separately, a sufficient number of valence quarks can always be 
added to the theory. 

4-4 Renormalization 

Equations (4.5)-(4.7) connect the spectral density to the basic field-theoretic correla¬ 
tion functions whose renormalization properties are well understood. The suggestion 
is then that the moments get renormalized through multiplication by the factor 
( Zp) n . Since the sum over the coordinates aq, ..., x n -\ in eq. (4.7) includes the 
short-distance regions, it is, however, not totally obvious that the moments thus 
renormalized will indeed have a well-defined continuum limit. 

In the continuum theory, and when inserted in correlation functions, the product 
Pi 2 (xi)P- 23 (x 2 ), for example, has a short-distance expansion 


(4.6) 


Pl2(xi)P23(x 2 ) 


r\j 


C(x 1 - X 2 )S 13(x 2 ) + • • • , 


(4.8) 



where the Wilson coefficient C(x) diverges like |a;| -3 (up to logarithms). This sin¬ 
gularity is integrable and does not give rise to any additional ultraviolet divergen¬ 
cies. As it turns out, all short-distance singularities of the correlation functions in 
eq. (4.7) are in fact integrable f. The ultraviolet divergencies of the moments are 
thus completely cancelled by the renormalization factor ( Zp) n and by the usual 
renormalization of the gauge coupling and the quark mass. 

Recalling the expansion (4.6) of the resolvent R(z) and its relation to the spectral 
density, it now follows that 


pB.( a ) — Zpp(Zpa) 


(4.9) 


has a universal continuum limit (once a particular renormalization condition for the 
pseudo-scalar density is adopted). The same must then also apply to the renorma¬ 
lized threshold 


<Tr — Z 


- 2 ; 


a, 


(4.10) 


since the threshold is a property of the spectral density. 

Without proof we mention in passing that the improved renormalized density in 
the 0(a)-improved theory is again given by eq. (4.9), provided Zp is replaced by 


Zp 


1 + bpam, q 
1 + bpparriq ’ 


(4.11) 


where 6p and 6pp are improvement coefficients and m q the additively renormalized 
bare quark mass (the notation is as in ref. [3]). The correction proportional to 6pp 
in this expression is in fact all that is needed to cancel the terms of order a that arise 
from the integrations over the short-distance singularities of the correlation function 
on the right of eq. (4.7). 

4-5 Relation to the current-quark mass 

The bare current-quark mass rri is usually determined through vacuum-to-pion ma¬ 
trix elements of the isovector axial current and density. From the renormalization 
properties of the matrix elements it then follows that the renormalized quark mass 


f When all coordinates x\,. ■ ■ ,x n —\ are scaled to zero, the operator product reduces to a coef¬ 
ficient function times the unit operator. The degree of divergence of the integral is 4 — n in this 
case, and convergence is thus guaranteed since n > 6. 
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is given by 


m R = Z A Z P 1 m, 


(4.12) 


where Z A denotes the axial-current renormalization constant. 

Once a definite renormalization condition for the isovector axial density is adopted, 
the threshold a R and the quark mass m R become physical quantities. In particular, 
in the continuum limit we have 

V® R _ ^ 
m R 

not only when the limit is reached from a lattice theory that preserves chiral sym¬ 
metry, because this ratio is dimensionless, unambiguously normalized and there¬ 
fore independent of the regularization. In terms of the bare quantities this implies 
= Z A rn. up to corrections of 0 (a) and independently of the normalization con¬ 
vention for the renormalized axial density [the renormalization constant Zp cancels 
in the ratio (4.13)]. 

As already mentioned, the spectral gap in finite volume may not be related in any 
simple way to the threshold of the spectral density in the thermodynamic limit. It is 
nevertheless instructive to compare the median fl of the gap distributions discussed 
in sect. 3 with the threshold Z A m. 


(4.13) 


Table 2. Comparison of the median /i with Zpjn 


Run 

Z A 

/ i/m 

/r — Z A m [MeV] 

(A) [MeV] 

Ai 

0.77(2)“ 

0.76(1) 

-0.5(5) 

2 . 8 ( 1 ) 

A-2 


0.75(1) 

-0.6(4) 

3.2(1) 

A 3 


0.80(3) 

0.6(5) 

3.5(1) 

a 4 


0.85(3) 

1.1(4) 

3.9(1) 

Bi 

0.78(2)“ 

0.77(1) 

-0.3(3) 

1.91(6) 

b 2 


0.79(1) 

-0.2(3) 

2.19(6) 

Ci 

0.77(2)“ 

0.72(1) 

-1.7(3) 

2.18(7) 


0.75(l) b 

0 . 68 ( 1 ) 

-5.4(5) 

2 . 2 ( 1 ) 


a RI-MOM method [22] 

b Schrodinger functional chiral Ward identity [23] 
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The numerically computed values of p/m quoted in table 2 actually agree quite 
well with the available estimates of Z&. There are, however, significant differences 
in the last two rows of the table (runs C\ and D i), which underlines the fact that 
there is currently no solid theoretical understanding of the gap distributions in finite 
volume. On the other hand, the absolute deviation of the median from the threshold 
is, in most cases, smaller than the average splitting (A) of the first four eigenvalues 
of |Qui | (see table 2; the figures quoted in the fourth column do not include the error 
on Za). In particular, the data are consistent with the working hypothesis on which 
our argumentation relied (cf. subsect. 4.1). 


5. Conclusions 

As explained in sect. 2, numerical simulations of the Wilson theory can be expected 
to be stable if the distribution of the spectral gap of the lattice Dirac operator is 
well separated from the origin. The range of stability may be defined through the 
inequality p > 3u, for example, where, as before, p and a denote the median and 
width of the distribution. Using the empirical relations p — Zm, and a — a/y/V, 
this bound becomes m > 3 a/ZyfV, which shows that the accessible range of quark 
masses depends on both the lattice spacing and the lattice size. 

Another form of the stability bound is obtained by multiplication with the ratio 
B = m^/2m, which is known to be practically independent of the quark mass [5,6]. 
For lattices of size 2 Lx L 3 , this leads to the inequality 

m n L > \J 3y/2aB / Z . (5.1) 

In particular, from our numerical studies of the (unimproved) Wilson theory, we 
deduce the stability bounds 

f 2.8 at a = 0.08 fm, 

m n L > < (5.2) 

l 2.3 at a = 0.06 fm, 

where the mass dependence of B and Z was neglected for simplicity. Note that the 
pion mass in these formulae is the one in infinite volume, and not the possibly much 
larger mass computed on a lattice of size L. 

It follows from these results that the range of stability includes all lattices where, 

say, a < 0.1 fm, L > 2 fm and m n L > 3. Simulations of the Wilson theory on such 
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lattices, using the known simulation algorithms, are thus expected to be safe from 
the problems mentioned in sect. 2. The 0(a)-improved theory is likely to behave in 
the same way, but so far this has only been checked on a single lattice with spacing 
a = 0.09 fm. If the median and width of the gap distribution are assumed to scale 
as in the unimproved theory, the stability bound deduced from the simulation of 
this lattice is m n L > 3.2. Extensive studies of the improved theory will evidently 
be required to confirm this result, which is very much in line with the bounds (5.2). 

The numerical simulations reported in this paper were performed on PC clusters 
at CERN, the Centro Enrico Fermi and the Institut fur Theoretische Physik der 
Universitat Bern (with a contribution from the Schweizerischer Nationalfonds). We 
are grateful to all these institutions for the continuous support given to this project. 
L. G. acknowledges partial support from the EU under contract HPRN-CT-2002- 
00311 (EURIDICE). 


Appendix A. Chebyshev accelerated subspace iteration 

On a finite lattice, the few lowest eigenvalues of A = Q ^ can be computed numeri¬ 
cally by minimizing the associated Ritz functional, for example [24,25]. This method 
is relatively tolerant of rounding errors, which is an important advantage on com¬ 
puters that do not support double-precision arithmetic. Otherwise there is a choice 
of algorithms that can be significantly faster. Subspace iteration with Chebyshev 
acceleration and eigenvector locking is one of them, and it is our aim, in the following 
paragraphs, to describe this method in some detail (see ref. [26], for example). 

A.l Power method 

A tight lower bound on the largest eigenvalue of A can be obtained by repeatedly 
applying the operator to a random quark field. Starting from some arbitrary field 
'(/’ with unit norm, the recursion 

X = M, -0 = x/llxll, (A.l) 

systematically enhances the upper spectral components of the field. The norm ||x|| 
then provides an increasingly accurate estimate of the largest eigenvalue. In practice 
20 iterations or so are usually sufficient for an accuracy better than 5%. 

When A is replaced by the shifted operator c — A, where c > ^||A|| is some fixed 
number, the power iteration converges to c minus the lowest eigenvalue of A. The 
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Fig. 7. Plot of the Chebyshev polynomial of order 12, showing the bounded oscil¬ 
latory behaviour in the interval [—1,1] and the rapid increase of the polynomial away 
from this range. 

latter can thus be calculated in this way. However, since the low eigenvalues of A 
are typically orders of magnitude smaller than c, accurate results are obtained only 
after a very large number of iterations. 

The power method is therefore not recommended for the computation of the low 
eigenvalues of A. unless it is combined with an acceleration technique. In particular, 
the shifted operator can be replaced by a Chebyshev polynomial and the iteration 
may be extended to a subspace of quark fields. Both of these modifications lead to 
significantly improved convergence rates. 

A.2 Chebyshev polynomials 

For \z\ < 1 the Chebyshev polynomials T 0 (z), Ti(z ),... are defined by 

Tk(z) = cos(kd), z = cos9. (A.2) 


They oscillate between —1 and +1 in this range and rapidly increase or decrease to 
±oo when \z\ > 1, depending on whether k is even or odd (see fig. 7). Through the 
linear transformation 


2 x — v — u 
v — u 


(A.3) 


the polynomials can be made to oscillate in an arbitrary interval x G [ u , u] instead 
of the standard range [—1,1]. 
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Fig. 8. Plot of the Chebyshev polynomial of order 50, scaled to the interval [u, v] = 
[0.1,49]. If the spectrum of A is contained in this range except for the few lowest 
eigenvalues (filled circles), the application of the polynomial to a given quark field 
enhances the field components along the subspace spanned by these modes. 

We may now replace x by A and choose [u, u] to contain the unwanted part of the 
spectrum of A. The other part is then strongly enhanced when the polynomial is 
applied to a given quark field (see fig. 8). Since the spectrum of the operator is not 
known beforehand, the interval bounds and the degree of the polynomial must be 
chosen adaptively in the course of the power iteration. 

A simple mathematical fact that will be used in this context is summarized by 

Lemma A.l. For any given even degree k, and real numbers x, v, 7 satisfying x < v 
and 7 > 1 , there exists a unique value of the interval bound u such that x < u < v 
and Tk(z) = 7 [where z is as in eq. (A.3)]. 

Proof: Since 7 is larger than 1, there is one and only one uj such that 

7 = coshcu, to > 0. (A.-4) 

The inequality x < u, on the other hand, implies z < — 1 and thus 

Tfc(z) = cosh(fczz), z = — coshzq (A.5) 

for some u > 0. We then conclude that the equation Tf,.(z) = 7 has a unique solution 
(respecting the specified constraints), which is given by u = oj/k. 
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Having computed z as a function of k and 7, the interval bound 


u = x + (v — 



(A.6) 


is obtained from eq. (A.3). This shows that u is uniquely calculable, and it is now 
also easy to verify that the expression (A.6) has all the required properties. 


A.3 Chebyshev accelerated subspace iteration 

The algorithm described in this subsection computes increasingly accurate approxi¬ 
mations to the n lowest eigenvalues of A and the associated eigenvectors. It operates 
on a d-dimensional subspace of quark fields, where d is usually taken to be quite a 
bit larger than n. 

The computation proceeds iteratively, starting from a random set of orthonormal 
quark fields ipi ,...., ipd- In each cycle of the iteration, the fields are updated one by 
one so that at the end of the cycle the whole basis of fields is replaced by a new one. 
More precisely, a cycle consists of the following four steps: 

1. Choose a Chebyshev polynomial T k (z) with even degree and an appropriate spec¬ 
tral range [u,v], taking the current estimates a\,... ,ctd of the eigenvalues as input. 

2. Update the basis vectors i = 1,..., d, one after another through 

m . . , 2 A — v — u , . 

X = T k (z)il>i, z = -, (A.7) 

v — u 

i— 1 

4> = X- (A-8) 

3 = 1 

A = 4>/\\4>\\- (A.9) 

The second and third equations here simply implement the modified Gram-Schmidt 
orthogonalization process. In particular, the new basis is guaranteed to be orthonor¬ 
mal. 

3. Rotate the fields among themselves so as to diagonalize the operator A in the 
subspace spanned by them, i.e. so that 

(V’i, Aipj) = SijOii for all i, j = 1 ,... d, 


Oi\ < Oi2 < • • • < Oid, 


after the transformation. 


19 


(A.10) 

(A.ll) 



4. Stop the algorithm if the approximate eigenvalues ai,... ,a n and the associated 
eigenvectors satisfy the chosen convergence criterion. 

This description is somewhat schematic and needs to be made more precise. It may 
not be obvious, for example, how to choose the polynomial in the first step, and 
there are various stopping criteria that may be applied. 


A.4 Choice of the Chehyshev polynomial 

Ideally the interval [u, v\ should contain the spectrum of A except for the d lowest 
eigenvalues. In particular, the upper limit v should be set to a value slightly larger 
than ||A||. If this number is not already known, it may be calculated at the beginning 
of the subspace iteration using the ordinary power method. 

To fix the lower bound u of the spectral interval and the degree k of the polynomial, 
a reasonable requirement is that 


T fcO)l*=ai =7 2 , T k (z)\ x=ad =7, (A.12) 

where 7 > 1 is the desired enhancement factor for the low modes (see fig. 8). In 
practice 7 = 3 appears to be a sensible choice, but trying other values of 7 may be 
worth while. 

Recalling lemma A.l, it is clear that the conditions (A.12) determine both u and 
k. Explicitly, if we define u and Co through 


coshu; = 7, coshu; = 

7 2 , 

(A.13) 

the conditions become 



u = ad + (v — ad) tanh 2 ( 

,2k)' 

(A.14) 

a d + {v - a d ) tanh 2 

= + (v — ai) tanh 2 ( 7 ^) ■ 

(A.15) 


Normally k is such that 2k 3> Co, and the expansion at large k of eq. (A. 15) then 
leads to 


k 


1 

2 


( v — 07) Co 1 — (v — ad) oj 2 
Old ~ Ol\ 


1/2 


(A.16) 


while u is determined by eq. (A. 14). Evidently k should be rounded to the closest 
even integer, and one may also wish to impose lower and upper limits on k at this 
point. 
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A.5 Stopping criterion 

An upper bound on the deviation of the calculated eigenvalues on from the exact 
eigenvalues of A may be obtained by computing the residues 

Pi = (A - aityi (i = 1,..., n) (A.17) 

and the maximal eigenvalue e 2 of the n x n residual matrix 

Rij = (pi,pj). (A.18) 

A well-known lemma then asserts that there are n orthonormal eigenvectors of A 
with eigenvalues on such that — &j| < e for all i = 1,..., n. 

This convergence criterion is safe but can be inefficient if the set {a\ ,... ,a n } of 
approximate eigenvalues divides into well separated subsets of one or more eigenval¬ 
ues. The error bounds obtained from the residual matrices associated to each subset 
of eigenvalues are then often quite a bit smaller than the bound obtained from the 
total residual matrix. Provided the subsets are indeed separated from one another, 
by a margin larger than the combined errors, these tighter bounds are completely 
safe too. 

A. 6 Eigenvector locking 

A fairly obvious property of the subspace iteration is that the lower eigenvalues con¬ 
verge faster than the higher ones. The algorithm can thus be accelerated somewhat, 
by locking the subsets of eigenvalues and eigenvectors that have already converged. 
Locking means that these eigenvectors are not updated in the second step of each 
subspace iteration cycle and that, in the third step, the operator A is diagonalized 
in the complementary subspace only. 

Another small acceleration is achieved by saving the last few eigenvectors, say 
ipd-r+ 1, • • •, V’dj to some auxiliary fields before they are updated in the second step. 
In the third step the saved fields may then be included in the Ritz diagonalization, 
i.e. A is diagonalized in a subspace of dimension d+r, but only the first d eigenvectors 
are kept after the diagonalization. 

A . 7 Rounding errors 

In order to avoid large rounding errors when the Chebyshev polynomials are applied 
to the quark fields, the Clenshaw recursion should be used [27]. Let ^ be a given 
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quark field and let us define 


Xj = 


2 A — v — a 

z = - 

v — u 


The computation then proceeds recursively according to 


(A.19) 


Xo =ip, Xi = (A.20) 

Xj+i Xj — 1 5 3 I? 2,... , (A.21) 

until the desired degree k is reached. 

The use of 32-bit arithmetic in the subspace iteration does not lead to uncontrolled 
rounding errors as long as the degree k of the Chebyshev polynomial is not too large. 
In general the significance loss in the Clenshaw recursion grows proportionally to 
k, and degrees below 100 or 200 may therefore be safe. On large lattices, however, 
the lowest eigenvalues of A tend to be closely spaced and many orders of magnitude 
smaller than the maximal eigenvalue. Polynomials with significantly larger degrees 
will be required under these conditions and the use of double-precision arithmetic 
then becomes indispensable. 
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